Review: kinds of spatial data

1) raster

Raster data are spatial data stored in a grid. Examples include Digital Elevation Models (DEM), spatially-interpolated climate or weather data (e.g., PRISM, DAYMET, WorldClim), or model predictions extrapolated across space (e.g., predictions of occupancy, abundance, density). The R package raster is all that is needed for most raster data handling in R.

2) vector

Vector data are data represented by points, lines, and polygons. Shapefiles, KMLs, and GeoJSONs are common formats used to represent geographic vector data. The R package sp handles most vector data needs in R. Other helpful functions can be found in the rgeos package, where tools similar to those found in the ArcGIS Analysis Toolbox (e.g., clip, union, buffer).

Plotting spatial data with Base R

First, we’ll load tmap, which we’ll use later to make some nice-looking maps. Here, we’re just going to load it because it has some “example” spatial data we can plot.

library(sp)
library(tmap)
data(World)
data(metro)

Both datasets, “World” and “metro,” are vector datasets. We’ll use Base R to plot these so you can see what they look like:

plot(World)

plot(metro, pch=16, cex=0.5)

Neat, but I want them on the same map. Let’s try this:

plot(World)
plot(metro, pch=16, cex=0.5, add=T)

I used the argument ‘pch=16’ to specify the symbol that I wanted to use for metro areas, a filled circle. More plotting characters can be found by googling it. The argument ‘cex=0.5’ controls the symbol size.

Our metropolitan areas didn’t show up. Why not? Usually it’s because the coordinate reference systems (CRS) of the datasets don’t match.

proj4string(World)
## [1] "+proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0"
proj4string(metro)
## [1] "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"

Herein lies a fundamental difference between making maps in R (especially base R) and making maps in ArcGIS. R doesn’t ‘project on the fly,’ so you’re responsible for lining everything up by getting all your datasets projected into a common CRS. The mapping functions in the tmap package have some project-on-the-fly functionality, but it’s best to get your layers in a common CRS anyway because most of the functions you will use to process spatial data demand it.

Reprojecting vector data

For this example, it doesn’t matter what CRS we use, just as long as all of our datasets share the same one. We’ll use the CRS from metro. To reproject vector data, we use the function sp::spTransform. The first argument is the dataset we want to reproject (“transform”), the second argument is the CRS, which needs to be cast into the CRS class. We’ll call this reprojected dataset “world” (lower-case w).

world <- spTransform(World, CRS(proj4string(metro)))

Now try plotting again:

plot(world)
plot(metro, pch = 16, cex = 0.5, add = T)

Base R is very flexible and you can make very sophisticated, visually attractive maps with it. It’ll just take a lot of code to get it looking right.

Raster data

Base R can also plot raster data. Let’s take a look at a common type of raster data that’s freely available: climate data. WorldClim has some basic data like average monthly precipitation and temperature, as well as a suite of “bioclimatic variables” that can be useful for predicting species distributions, etc. We’ll download the bioclimatic variables and look at a few of them.

library(raster)
paths <- list.files("~/Downloads/wc2.0_10m_bio", ".tif", full.names = T)
raster.list <- lapply(paths, function(x) raster(x))
bio <- do.call(brick, raster.list)
## Warning in as.POSIXlt.POSIXct(x, tz): unknown timezone 'zone/tz/2018c.1.0/
## zoneinfo/America/Denver'
class(bio)
## [1] "RasterBrick"
## attr(,"package")
## [1] "raster"

We now have a “RasterBrick” containing all the bioclimatic variables as layers or bands.

nlayers(bio)
## [1] 19
names(bio)
##  [1] "wc2.0_bio_10m_01" "wc2.0_bio_10m_02" "wc2.0_bio_10m_03"
##  [4] "wc2.0_bio_10m_04" "wc2.0_bio_10m_05" "wc2.0_bio_10m_06"
##  [7] "wc2.0_bio_10m_07" "wc2.0_bio_10m_08" "wc2.0_bio_10m_09"
## [10] "wc2.0_bio_10m_10" "wc2.0_bio_10m_11" "wc2.0_bio_10m_12"
## [13] "wc2.0_bio_10m_13" "wc2.0_bio_10m_14" "wc2.0_bio_10m_15"
## [16] "wc2.0_bio_10m_16" "wc2.0_bio_10m_17" "wc2.0_bio_10m_18"
## [19] "wc2.0_bio_10m_19"

Maybe you’re interested in precipitation. According to the WorldClim website:

  • BI012 = Annual Precipitation
  • BIO13 = Precipitation of Wettest Month
  • BIO14 = Precipitation of Driest Month
  • BIO15 = Precipitation Seasonality (Coefficient of Variation).

These correspond to layers 12 through 15, so let’s subset our RasterBrick to just those layers and take a look at them.

precip <- bio[[12:15]]
class(precip)
## [1] "RasterStack"
## attr(,"package")
## [1] "raster"
plot(precip)

The colors are alright, but they can be easily changed. The function colorRampPalette is a convenient and flexible way to make interpolated color ramps for continuous data like these.

precip.colors <- colorRampPalette(c("#fc8d59", "#ffffbf", "#4575b4"), bias = 4)
plot(precip, col = precip.colors(100))

Reprojecting raster data in R

Reprojecting raster data requires a different function than that used on vector data. Unsurprisingly, it’s found in the raster package, raster::projectRaster. The same function is used whether you’re dealing with a single layer/band (a RasterLayer) or a multi-layer RasterStack or RasterBrick.

precip <- projectRaster(precip, crs = CRS(proj4string(world)))

Now let’s add the country boundaries and metro areas,

plot(precip[[1]], col = precip.colors(100))
lines(world, lwd = 0.25, col = "grey50")
points(metro, cex = 0.8 * (metro$pop2010/max(metro$pop2010)))

Getting serious: making better maps with tmap

tmap works much like ggplot.

tm_shape(precip) +
  tm_raster(col="wc2.0_bio_10m_12", title="Annual precipitation", palette=precip.colors(100), n=100, legend.show = F) +
tm_shape(world) +
  tm_borders("grey35", lwd=0.5) +
tm_shape(metro) +
  tm_symbols(size="pop2010", title.size="Metro population", col="black", shape=1, scale=0.5)

Now if we want to plot the same basic map but display multiple attributes in different panels, just pass a vector of raster layer names to the “col” argument in tm_raster().

Recall that these can be accessed with names():

names(precip)
## [1] "wc2.0_bio_10m_12" "wc2.0_bio_10m_13" "wc2.0_bio_10m_14"
## [4] "wc2.0_bio_10m_15"
tm_shape(precip) +
  tm_raster(col=names(precip), palette=precip.colors(100), n=100, legend.show = F) +
tm_shape(world) +
  tm_borders("grey35", lwd=0.5) +
tm_shape(metro) +
  tm_symbols(size="pop2010", title.size="Metro population", col="black", shape=1, scale=0.5)

To change aspects of the general appearance of the map, use the function tm_layout

tm_shape(precip) +
  tm_raster(col=names(precip), palette=precip.colors(100), n=100, legend.show = F) +
  tm_facets(free.scales=T) +
tm_shape(world) +
  tm_borders("grey35", lwd=0.5) +
tm_shape(metro) +
  tm_symbols(size="pop2010", title.size="Metro population", col="black", shape=1, scale=0.5) +
tm_layout(bg.color="grey80",
          legend.bg.color="white",
          legend.outside=T,
          legend.outside.position="right",
          panel.labels=c("Annual Precipitation",
                         "Precipitation of Wettest Month",
                         "Precipitation of Driest Month",
                         "Precipitation Seasonality (CV)"),
          panel.label.bg.color="white")

Say we want to change the frame of our map to only show Africa. This can easily be done by specifying the bbox (bounding box) using tmaptools::bb(). First, let’s have a look at the world dataset to see how we can subset this SpatialPolygonsDataFrame to Africa.

head(world)
##    iso_a3                 name           sovereignt     continent
## 2     AFG          Afghanistan          Afghanistan          Asia
## 3     AGO               Angola               Angola        Africa
## 5     ALB              Albania              Albania        Europe
## 8     ARE United Arab Emirates United Arab Emirates          Asia
## 9     ARG            Argentina            Argentina South America
## 10    ARM              Armenia              Armenia          Asia
##          subregion    area  pop_est pop_est_dens gdp_md_est gdp_cap_est
## 2    Southern Asia  652860 28400000     43.50090      22270    784.1549
## 3    Middle Africa 1246700 12799293     10.26654     110300   8617.6635
## 5  Southern Europe   27400  3639453    132.82675      21810   5992.6588
## 8     Western Asia   83600  4798491     57.39822     184300  38407.9078
## 9    South America 2736690 40913584     14.95003     573900  14027.1261
## 10    Western Asia   28470  2967004    104.21510      18770   6326.2469
##                      economy              income_grp life_exp well_being
## 2  7. Least developed region           5. Low income     48.7   4.758381
## 3  7. Least developed region  3. Upper middle income     51.1   4.206092
## 5       6. Developing region  4. Lower middle income     76.9   5.268937
## 8       6. Developing region 2. High income: nonOECD     76.5   7.196803
## 9    5. Emerging region: G20  3. Upper middle income     75.9   6.441067
## 10      6. Developing region  4. Lower middle income     74.2   4.367811
##         HPI
## 2  36.75366
## 3  33.20143
## 5  54.05118
## 8  31.77827
## 9  54.05504
## 10 46.00319

Conveniently, there’s a column called “continent” to which we can apply a logical statement to subset the data, like so:

tm_shape(precip, bbox=tmaptools::bb(world[world$continent == "Africa",], ext=1.1)) +
  tm_raster(col=names(precip), palette=precip.colors(100), n=100, legend.show = F) +
  tm_facets(free.scales=T) +
tm_shape(world) +
  tm_borders("grey35", lwd=0.5) +
tm_shape(metro) +
  tm_symbols(size="pop2010", title.size="Metro population", col="black", shape=1, scale=0.75) +
tm_layout(bg.color="grey80",
          legend.bg.color="white",
          legend.outside=T,
          legend.outside.position="right",
          panel.labels=c("Annual Precipitation",
                         "Precipitation of Wettest Month",
                         "Precipitation of Driest Month",
                         "Precipitation Seasonality (CV)"),
          panel.label.bg.color="white")

Getting interactive: tmap view mode

Note: the following is modified only slightly from the vignette “tmap-modes” (see vignette(“tmap-modes”, package=“tmap”))

Static maps are fine, but the world is going interactive (for example, see this this or this).

There are a growing number of tools that facilitate the production of interactive figures, including maps, in R. The tmap package uses Leaflet to implement interactive maps; this feature is accessed by switching into “view mode.”

First we’ll make a static map, but this time we’ll assign it to an object instead of simply printing it.

metro$growth <- (metro$pop2020 - metro$pop2010) / (metro$pop2010 * 10) * 100

mapWorld <- tm_shape(world) +
    tm_polygons("income_grp", palette="-Blues", contrast=.7, id="name", title="Income group") +
    tm_shape(metro) +
    tm_bubbles("pop2010", col = "growth", 
               border.col = "black", border.alpha = .5, 
               style="fixed", breaks=c(-Inf, seq(0, 6, by=2), Inf),
               palette="-RdYlBu", contrast=1, 
               title.size="Metro population", 
               title.col="Growth rate (%)", id="name") + 
    tm_style_gray() +
    tm_format_World()

To plot, just print the declared variable mapWorld:

mapWorld

The current mode can be obtained and set with the function tmap_mode:

tmap_mode("view")
## tmap mode set to interactive viewing

Now the mode is set to view, we can interactively view it by printing it again:

mapWorld
## Legend for symbol sizes not available in view mode.

In interactive mode, a leaflet widget is created and shown. If you want to change or extend the widget, you can use the function tmap_leaflet to obtain the leaflet widget object, and use leaflet’s own functions to adjust it.